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ABSTRACT 



Aims. We measure the turbulent resistivity in the nonlinear regime of the MRI, and evaluate the turbulent magnetic Prandtl number. 
Methods. We perform a set of numerical simulations with the Eulerian finite volume codes Athena and Ramses in the framework of 
the shearing box model. We consider models including explicit dissipation coefficients and magnetic field topologies such that the net 
magnetic flux threading the box in both the vertical and azimuthal directions vanishes. 

Results. We first demonstrate good agreement between the two codes by comparing the properties of the turbulent states in simulations 
having identical microscopic diffusion coefficients (viscosity and resistivity). We find the properties of the turbulence do not change 
when the box size is increased in the radial direction, provided it is elongated in the azimuthal direction. To measure the turbulent 
resistivity in the disk, we impose a fixed electromotive force on the flow and measure the amplitude of the saturated magnetic field 
that results. We obtain a turbulent resistivity that is in rough agreement with mean field theories like the Second Order Smoothing 
Approximation. The numerical value translates into a turbulent magnetic Prandtl number Pm, of order unity. Pm, appears to be an 
increasing function of the forcing we impose. It also becomes smaller as the box size is increased in the radial direction, in good 
agreement with previous results obtained in very large boxes. 

Conclusions. Our results are in general agreement with other recently published papers studying the same problem but using different 
methodology. Thus, our conclusion that Pm, is of order unity appears robust. 

Key words. Accretion, accretion discs - MHD - Methods: numerical 



1. Introduction 

Magnetic fields play a central role in accretion disk 
theory. Thro ugh the magn etorotational instability (MRI, 
iBalbus & Ha wlev 19911 1 19981) . they can explain why the flow 
is turbulent and how angular momentum is transported outward. 
They are als o thought to be responsible for jet launching and 
coUimation dBlandford & Pavndl 19821 iPudritz et al]|2007l) . For 
both phenomena, whether or not a mean magnetic flux threads 
the disk has important consequences. Numerical simulations 
have indeed shown that the saturation level of the MRI depends 
strongly on the net flux: the famous a parameter (a measure of 
the rate of angular mome ntum transport through the disk, see 
IShakura & Sunvaevll 19731) is directly proportional to the square 
of the vertical magnetic field dHawlev et alJll995l) . Similarly, a 
strong vertical field is mandator y for accretion-ejectio n models 
to be efficient at launching jets ( Casse & Ferreirall200 0l). 

However, MHD turbulence not only can transport angu- 
lar momentum outward in the disk, but also can diffuse mag- 
netic flux radially. This possibility questions the very presence 
of a mean magnetic flux in the inner part of accretion disks 
dvan Ballegooiienl Il989t: iLubow et all 1 19941) . Ultimately, how 
much magnetic flux resides in the inner disk depends on a com- 
petition between the effectiveness of outward angular momen- 
tum transport (i.e. how well mass and magnetic flux are advected 
inward) and magnetic flux diffusion. The former effect can be 
identified with a "turbulent viscosity". Similarly, the importance 
of the latter can be assessed through an equivalent "turbulent re- 
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sistivity" rjt which has physical consequences similar to that of 
a microscopic resistivity. The purpose of this paper is to mea- 
sure rj, accurately by using a set of local numerical simulations 
in which the flow is turbulent because of the MRI. 

We note that tw o recently published papers have considered 
the same problem dGuan & Gammid 120091; iLesur & Longarettil 
2009). Although all such studies (including our own) have used 
numerical simulations in the local shearing box model, there 
are significant differences between the approach used in each 
case. First, the numerical methods used to sol ve the MHD equa- 
tions are different in all three approaches: iGuan & Gammid 
d2009h use d numerical methods similar to those in the 
ZEUS code dSt one & Norman 1992; Hawlev & Stonell995h and 
iLesur & Longaret ti (2009) used a pseudo-spectral code in the in- 
compres sible limit. By contrast, we use two finite volume co des, 
Athena ("Gardiner & Stone"2005a', '2008"; 'Stone etalj|2008') and 
Ramses (Teyssier 2002; Fromang et al. 2006). In addition, the 
method used in this paper to measure the turbulent resistivity 

i see se ctio n ^ is different froin that u sed by Guan & Gammid 
2009h and ILesur & Longarettil d2009l) . The former add a large 
scale magnetic field to an already turbulent flow and associate 
the rate of decay of this field to the turbulent resistivity in the 
disk. The latter impose at all times an additional magnetic field 
in the computational domain, measure the time averaged elec- 
tromotive force (EMF) associated with this field, and use the 
amplitude of the EMF to measure the turbulent resistivity. As 
described below, our approach is to impose a forcing EMF on 
the flow and relate the time averaged structure of the saturated 
field that results to the value rj,. Another difference between our 
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work and previous studies is that we study turbulence driven by 
a magnetic field with not net flux threading the box in any di- 
rection, while Guan & Gammie (2009) and Lesur & Longare^tl 
(l2009l) perform numerical simulations in the presence of either 
a net vertical or a net toroidal field. It is well known that the 
presence of a net field affects the properties of the turbulence 
that develops (iHawley et al.lll995t iGuan et al.1 12009 ). Thus, an 
additional goal of this paper will be to study how much the final 
results obtained in all these studies depends on the details of the 
simulation setup. This will be useful to assess the robustness of 
the results. 

The plan of this paper is as follows. In section |2] we detail 
our method to measure the turbulent resistivity and the numer- 
ical setup we used. Section |3] serves essentially to validate our 
numerical codes by considering the case of a vanishing forc- 
ing EMF (in which case the simulations are standard shear- 
ing box simulations). We show that we obtain turbulent flows 
whose time- and volume-averaged properties are identical with 
both Athena and Ramses. We also study whether these results 
are modified when we consider boxes extended in the r adial di- 
recti on, as suggested recently in the literature (Johanse n et al.l 
12008 ). These simulations serve as a basis for section|4]in which 
we measure the turbulent resistivity in a number of situations, 
varying the strength and direction of the forcing EMFs, and 
the size of the box. Finally, in section |5] we summarize our 
res ults, compare them with those of Guan & Gammie (2009) 
and lLesur & Longarettil (l2009h . and stress the limitations of our 
study. 



2. The setup 

2.1. Equations 

We solve the equations of magnetohy drodynamics (MHD) in the 
frame work of the shearing box model (Goldreich & Lvnd en-Bell 
1965|) using two different cod es: Athena (Gardiner & Stone 
^2005a, 2008; Sto ne et all l2008l) and Ramses dTevssieij 120021 
Fromang et al.. Explicit microscopic dissipation, namely 

viscosity v and resistivity rj, are included in the simulations 
(note the microscopic resistivity rj should not be confused with 
the turbulent resistivity i], that will be introduced below). It 
has been shown in the last few years that both are impor- 
tant i n setting the satu rated stat e of turbulence dri ven by the 
MRI dLesur & Longare tti 2007; F romang et alj2007h . We use an 
isothermal equation of state P = pc^ to relate density p and pres- 
sure P, where cq is the sound speed. In a Cartesian coordinate 
system {x,y,z) with unit vectors (i, j, k) in the radial, azimuthal 
and vertical directions respectively, the equations for mass and 
momentum conservation are: 



^ + V.(pv) = 

of 



(1) 



dpv 

— + V-(pvv - BB) + VP,o, = 3pQxi - 2p£l xv + V-T, (2) 

ot 

where Q. is the angular velocity of the shearing box around 
the central object, B is the magnetic field, Ptoi = P + 
(B-B) /2 is the total pressure, and T is the viscous stress ten- 
sor dLandau & Lifshitj Il959i) . The first two source terms on 
the right hand side of Eq. (|2]l represent the tidal gravity and 
Coriolis forces. Their implementation in Athena and Ramses 
is n ot straightforward. We followed the method described by 
[Gardiner & Stone (2005b) (see also Stone & Gardiner (2009)) 
and use a Crank-Nicholson algorithm to update the momentum 
fluctuations in time due to these source terms. The divergence of 



the viscous stress tensor in Eq. (|2| is differenced in the conser- 
vative form so as to exactly conserve total momentum. 

In order to study turbulent resistivity, we use a modified form 
of the induction equation: 



dB 

— = Vx (vxB - nVxB + E) 

at 



(3) 



In this form we have introduced E, an imposed electromotive 
force (EMF), in addition to the usual inductive and resistive 
terms. When £ = 0, as in section|3] we recover the usual form of 
the induction equation for resistive MHD. In section|4] we con- 
sider solutions to Eqs. ([TJ-([3]l using the following expression for 
the imposed EMF: 



(2nx \ I 2nz 

— Ij + ^o.i cosl — 



(4) 



where Ae^x and Ae^, arc the radial and vertical wavelengths for 
the forcing EMF, respectively. To gain insight into the effect of 
E, it is instructive to consider the case v = 0. In this situation, 
the magnetic field reaches a steady state in which forcing due to 
E is balanced by resistive dissipation, with the field given by the 
following solution 



B 



2nri 



2nz\ . 



EOyAE,X 

2nri 



'2nx\ 
sin I k . 



(5) 



Eq. ^ shows that the effect of the radial (x) component of the 
EMF is to create a vertically varying azimuthal magnetic field, 
while the effect of the azimuthal (y) component of the EMF is to 
create a radially varying vertical magnetic field. The amplitude 
of the field is determined by the amplitudes and wavelengths of 
the forcing and the resistivity. In a turbulent flow we expect the 
same balance will be achieved: turbulent velocity fluctuations 
will diffuse and reconnect the magnetic field that the driving 
EMF is trying to build. By analogy with Eq. (|5]l, in a turbulent 
flow, we can define a turbulent resistivity t], by measuring the 
amplitude of the spatially varying field resulting from the forc- 
ing. More specifically, we will consider two cases. First, Eq^ - 
and £■() V > 0, in which case the forcing results in a purely vertical 
field of amplitude Z?". This allows us to measure the turbulent re- 
sistivity rii x in the radial direction through the following relation, 
obtained using Eq. (|5]l. 



Eoy Ae,x 

2^'s9' 



(6) 



Similarly, the case Eqx > and Eqx = yields a purely azimuthal 
field varying with z and with an amplitude B° that measures the 
turbulent resistivity 77, ^ in the vertical direction according to 



Eqx '^E.z 

2^'so' 



(7) 



2.2. Numerical metiiod 

We solved Eq. ([U, (|2|i and (O using Athena and Ramses in a 
box of size (Lx,Ly,L,) = {H,4H,H), where H - cq/Q defines 
the vertical scale height (thickness) of the disk. We also used 
Athena to compute a few runs with a larger extent in the radial 
direction: (Lx, L, , L~) = (4//, 4H, H). It has been shown recently 
by Johnson et al. (2008) that truncation error in the shearing box 
can create numerical artifacts, such as a density minima at the 
center of the box (where the shear velocity with respect to the 
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grid is zero). We have implemented an orbital advection algo- 
rithm ('Masset"2000'; 'Gammie"2001'; 'Johnson & Gammi^|2oio5l 
[Johnson et al. 2008) in Athena (Stone & Gardiner 2009) to alle- 
viate this problem, and at the same time to allow for much larger 
timesteps in wide radial boxes. 

All of our simulations are performed with no net magnetic 
flux in either the vertical or azimuthal direction. At f = 0, we 
start with a purely vertical magnetic field varying sinusoidally 
with X so that the mean field vanishes. Random velocity pertur- 
bations of small amplitude are applied to the background state, 
consisting of a uniform density gas with a linear shear profile: 
V = (0,-3Qjc/2,0). We adopt Q = cq = 10^1 In the rest 
of this paper, we measure time in units of the orbital period 
Torb - 2n/Q.. In all of the simulations we present below, we 
set V and rj such that the Reynolds number Re = cqH/t] = 3125 
and the magnetic Prandtl number Pm = 4. Th ese coefficients 
are ide ntical to the run labeled 128Re3125Pm4 of lFromang et alj 
(l2007 l). and are known to lead to sustained turbulence over large 
numbers of orbital times. For explicit dissipation to dominate 
over numerical dissipation, we used a resolution of 128 grid 
points per H, or {Nx,Ny,Nz) - (128,512, 128) for calculations 
that span one scale height in the radial direction. Such a res- 
olution has been shown to give resolved solutions when using 
2nd order Eulerian codes like ZEUS (Fromang et al. 2007) for 
micr oscopic dissip ation coefficients of the magnitude adopted 
here. ISimon et al.l (,2009) have recently shown that the numeri- 
cal dissipation in Athena, likely comparable to Ramses, is sim- 
ilar in magnitude to that of ZEUS when performing numerical 
simulations of MRJ-induced MHD turbulence in the shearing 
box. Thus, we expect the same resolution needed with ZEUS is 
also appropriate for the runs presented here that use Athena or 
Ramses. In simulations with larger boxes, we scaled the resolu- 
tion accordingly, (Nx,Ny,Nz) - (512,512, 128), so that the cell 
size remains the same. 

Regardless of the code we use or the size of the computa- 
tional domain, to study the turbulent resistivity induced by the 
MRI our strategy is as follows: we first performed a run without 
forcing (i.e. £ = 0) to let the MRI develop and for turbulence 
to reach a steady state. Typically ffiese runs are evolved for 100 
orbits. They also serve as a useful comparison between codes, 
and to compare small vs. large boxes. At f = 30, we restarted 
the simulations with forcing of various amplitudes imposed, as 
described above. With forcing, the flow evolves toward a new 
quasi steady-state that we use to evaluate the turbulent resistiv- 
ity induced by the MRI. Before moving on to a description of the 
results with forcing in section|4] we first describe the properties 
of MHD turbulence as obtained in runs without forcing in the 
next section. 



3. Runs without forcing 

The parameters of the runs we present in this section, along with 
time averaged values from selected quantities measured from 
the runs, are given in Table [1] Each run is labeled according 
to the code used ("R" for Ramses, or "A" for Athena) and the 
size of the box: "SB" for Small Box, i.e. simulations performed 
wiffi (Lx,Ly,Lz) - (H,4H,H) and "LB" for Large Box, i.e. 
(Lx,Ly,Lz) = (4H,4H, H). The resolution and duration of the 
runs are given in columns 3 and 4, respectively. In column 5, 
we report for each model the value of a time averaged between 
f = 50 and the end of the simulation. As is usual, it is defined 
by the sum of the Reynolds and Maxwell stress tensors, T^^^ 



and T^"^, normaUzed by the volume averaged thermal pressure 
<P>: 

a - - '■ — , (8) 

<P> <P> ^ ' 

where v~v and v\ are y and z averages of Vj and Vy respectively, 
and <> denotes a volume average. Finally, the last three columns 
in Table [1] give the amplitude of the gas velocity fluctuations 
along the three spatial coordinates. 

We performed three runs. Two of them, R-SB and A-SB, 
share identical parameters but use different codes, Ramses and 
Athena respectively, in order to compare the results from both 
codes. The third one, A-LB, was performed with Athena in a 
larger domain. The last four columns of Table [T] show that the 
saturated state of the turbulence is almost identical in the three 
models. For example, a = 0.014, 0.015 and 0.014 respectively 
in model R-SB, A-SB and A-LB. This agreement is confirmed 
by the three panels of figure [T] which correspond, from left to 
right, to R-SB, A-SB and A-LB. All show the time history of 
r^'V < P > (lower dashed line), T^""/ < P > (upper dashed 
line) and a (solid line). The curves obtained in the three sim- 
ulations are very similar and confirm the time averaged a val- 
ues given in Table [1] These results have two implications. First, 
the good agreement between R-SB and A-SB gives confidence 
in both codes, as the two models have identical parameters and 
differ only in the algorithms. Thus, it validates both the imple- 
mentation of the shearing box model in Ramses and the imple- 
mentation of orbital advection in Athena. Given the rms fluctu- 
ations in a reported in Table [1] these results are also consistent 
with the value of a ~ 0.01 quoted by iFromang et al.l (l2007l) for 
their model 128Re3125Pm4. We note, nonetheless, that ffie re- 
sults obtained here with Athena and Ramses appear to lead to 
slightly larger value of the angular momentum transport. It is 
unclear (and beyond the scope of this paper) whether this dif- 
ference is significant. It may be partly due to our box size being 
slightly larger i n the azimuthal direction compared to that used in 
Froman g et al.l (12007). or it could be due to small differences in 
the numerical dissipation between the various codes. The second 
result that emerges from the models presented in this section is 
the insensitivity of the turbulence properties to the box size. The 
time averaged value of a in model A-LB is almost identical to 
the other two. The time history presented on the right panel of 
figure [1] also shows good agreement with the other two panels. 
The only difference is that the fluctuations around the mean a 
value are smaller in the case of the big box. This translates into 
a standard deviation which is about three times smaller in A-LB 
than in the other two models. Again, the reasons for this differ- 
ence are unclear. It may be due to better statistics in the large 
box (simply due to the larger volume of the simulations) that 
average over extreme events. Alternatively, it might be due to 
ffie lar ger number of parasitic mo des than can develop in larger 
boxes jPessah & Goodmanl2009h or to a more efficient nonlinear 
coupling betwe en the more nume rous turbulent modes present 
in a larger box dLatter et al.ll2009l) . Both effects would reduce 
the lifetime (and therefore the influence) of channel modes. The 
similarity between the small and large boxes model is further 
confirmed by figure |2] and [3] The former shows snapshots of the 
density (left panel). By (middle panel) and (right panel) in 
the (x, z) plane for model R-SB at f = 70 (very similar figures 
are obtained for model A-SB). The equivalent figures for model 
A-LB at f = 70 are shown in figure [3] The figure shows that 
the structure of the flow in the large box is very similar to four 
patches of the small box model repeated next to each other. 
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Model 


Box size 


Resolution 


tlim 




a 


<Jv,./co 


(5v,,/co 


5v-Jcq 


R~SB 
A-SB 
A-LB 


(//, 4H, H) 
(H, AH, H) 
{AH, AH, H) 


(128,512, 128) 
(128,512, 128) 
(512,512, 128) 


160 
100 
100 


1.4 X 10- 

1.5 X 10- 
1.4 X 10- 


± 3.6 X 10--' 
-2 ± 5.7 X 10-3 
-2 + 1.4x 10-3 


0.0982 
0.0936 
0.0945 


0.0776 
0.0774 
0.0815 


0.0609 
0.0588 
0.0573 



Table 1. Properties of MHD turbulence when forcing is turned off. The first column labels the model. The box size, resolution and 
duration (in orbits) of the simulations are given in columns two, three and four The time averaged properties of the turbulence, i.e. 
the value of a (total stress normalized by the thermal pressure) and the velocity fluctuations along the three spatial (normaUzed by 
the speed of sound) appear in the following four columns. 




Fig. 2. Structure of the flow in the (x, z) plane obtained in model R-SB at t=70 orbits. The left, middle and right panels show the 
density, azimuthal component of the magnetic field, and vertical component of the velocity respectively. Similar results are obtained 
in model A-SB. 



Other features typical of shearing box simulations of 
the MRI, Uk e density wave propagating radially in the box 
jHeinemann & Papaloizou 2008a b), are clearly seen in the den- 
sity field regardless of the box size. Finally, as expected for a 
simulation having Pm > 1, the smallest scale structure seen in 
the velocity field is generally larger than the smallest scale struc- 
ture seen in the magnetic field. The latter shows elongated fila- 
ments reminisce nt of high Pm small scale dynamo simulations 
(see for example [Schekochihin et alj[2007l) . 

Recently, iJohansen et alj (l2008h also reported numerical 
simulations of MHD turbulence in the shearing box in large do- 
mains. They found an increase of a when going from small (//) 
to large (AH) boxes. Their result differs from those reported here, 
where we find no sensitivity of a on box size. The reason proba- 
bly lies in the fact that we always use domains with an azimuthal 
extend of AH, even when th e radial domain is sma ll (i.e. when 
Lx - H). On the other hand. iJohansen et all (l2008l) used square 
domains, with - Ly - H in their small box. It is well known 



that the azimuthal extent of the domain affec ts the saturated level 
of the turbulence dHawley et al.lll995L[T99^ . We conclude that, 
for fixed and large azimuthal extent of the domain, the converged 
value of a is insensitive to the radial dimension of the box, at 
least for the values of the Reynolds and Prandtl numbers studied 
here. 

Another property of the flow reported by IJohansen et alJ 
(2008.) is the existence of large scale, long-lived zonal flows that 
generate axisymmetric density features in large boxes. We have 
looked for, and found, such features in our large box simula- 
tions as well. An example is shown in figure |4] a space-time 
dia gram of the density i n model A-LB. It is similar to figure 6 
of J ohansen et al.l (l2008i) and shows a similar amplitude and life- 
time for large scale density features, indicating that a zonal flows 
also develops in our simulations. As both results were obtained 
using complet ely different codes, our results confirm the conclu- 
sion drawn by IJohansen et al.l ( l2008l) : zonal flows appear to be a 
robust feature in MRI-driven turbulence in the shearing box. 
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Fig. 4. Space-time diagram showing the radial variation of the density, averaged along y and z, in model A-L B. It shows the 
appear ance of long-lived density features. They are associated with large scale zonal flows as recently reported bv lJohansen et al.l 
(I2008h . 



This completes our description of the models computed in 
the absence of a forcing EMF. The flow in these models pro- 
vides a starting point to study turbulent resistivity induced by the 
MRI. Namely, at f = 30, we restarted the models using various 
amplitudes and directions for the forcing EMF, E. The following 
section describe these results. 



4. Turbulent resistivity 

The properties and results of the runs we performed are listed 
in Table |2] The first column gives a label for the name, coded 
according to the following convention: "direction-amplitude- 



code-parameters" where "direction" can be either "Ey" or "Ex" 
and indicates the direction of the forcing EMF (we considered 
only forcing aligned along either the radial or the azimuthal di- 
rection, respectively creating vertical or toroidal field), "ampli- 
tude" gives the amplitude of E, "code" reports the code we used 
for that model ("R" for Ramses, and "A" for Athena), while "pa- 
rameters", whenever present, gives additional parameters that 
will be specified when needed in the text. For example, model 
Ey-4E-10-R was performed with Ramses using a forcing EMF 
in the azimuthal direction of amplitude 4 x 10 The other pa- 
rameters reported in Table|2]are the duration over which the data 
were averaged (column 2), the amplitude of the forcing EMF 
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Model 




\E\ 




a 




SvJcq 




Sv-Jca 




1i 




Pm, 


Fv-4F- 1 0-R 


50 


4 X 10" 


10 


2.3 X 10" 


-1 


0.1214 


0.1060 


0.0789 


5.6 X \(f 


1.07 X 10" 


-5 


1.47 


Ey-4E-10-R-LONG 


120 


4 X 10" 


10 


2.0 X 10" 


-2 


0.1173 


0.1001 


0.0758 


5.5 X 10* 


1.06 X 10" 


-5 


1.27 


Ey-lE-9-R 


50 


1 X 10 


-9 


2.6 X 10" 


-2 


0.1344 


0.1231 


0.0896 


1.7 X W 


1.46 X 10" 


-5 


1.21 


Ey-4E-9-R 


50 


4 X 10 


-9 


5.2 X 10" 


-2 


0.1831 


0.1885 


0.1246 


1.6 X 10' 


1.81 X 10" 


-5 


1.91 


Ey-8E-9-R 


50 


8 X 10 


-9 


7.6 X 10" 


-2 


0.2209 


0.2458 


0.1493 


7.3 X 10^ 


2.43 X 10" 


-5 


2.10 


Ey-4E-10-A 


50 


4 X 10" 


10 


1.7 X 10" 


-1 


0.1020 


0.0846 


0633 


5.8 X lO'' 


1.08 X 10" 


-5 


1.02 


Ey-4E-9-A 


50 


4x 10 


-9 


5.3 X 10" 


-2 


0.1766 


0.1876 


0.1189 


1.2 X 10' 


1.56 X 10" 


-5 


2.28 


Ey-4E-10-A-L4-lx4 


48 


4x 10" 


10 


2.6 X 10" 


-1 


0.1242 


0.1202 


0.0798 


2.8 X lO'' 


3.02 X 10" 


-5 


0.57 


Ey-4E-9-A-L4-lx4 


50 


4x 10 


-9 


1.0 X 10" 


-1 


0.2371 


0.2624 


0.1555 


1.4 X 10' 


6.68 X 10" 


-5 


1.00 


Ey-4E-9-A-L4-lxl 


50 


4x 10 


-9 


5.3 X 10" 


-2 


0.1789 


0.1875 


0.1204 


1.6 X 10' 


1.82 X 10" 


-5 


1.93 


EX-4E-9-R 


50 


4x 10 


-9 


2.4 X 10" 


-2 


0.1292 


0.1162 


0.08476 


5.3 X 10- 


1.05 X 10" 


-5 


1.57 


EX-8E-9-R 


50 


8x 10 


-9 


4.6 X 10" 


-2 


0.1746 


0.1747 


0.1160 


1.3 X 10^ 


1.01 X 10" 


-5 


3.00 



Table 2. Properties of MHD turbulence when a forcing EMF is turned on. The first column gives the model label, the second the 
duration over which the data are averaged, and the third the amplitude of the forcing EMF used in that case. Subsequent columns 
give time averaged values of various quantities of physical interest, namely the mean values of a, the maximum value of (3 (ratio 
between the thermal and magnetic pressure), the value of the turbulent resistivity that results from fitting the mean vertical magnetic 
field, and the turbulent magnetic Prandtl number of the flow. 



(column 3), the value of a (column 4) and the amplitude of the 
velocity fluctuations along the three spatial coordinates (column 
5, 6 and 7). 

The forcing EMF we imposed on the flow results in a steady- 
state magnetic field profile. We fitted the profile (averaged in 
time over the duration of the run and in space over the two direc- 
tions perpendicular to the direction of variation) by a sinusoidal 
function. The amplitude of the fit, B fit, is expressed via the pa- 
rameter 



<p> 



(9) 



and given in column 8 on Table |2] The turbulent resistivity rj, 
responsible for that amplitude is calculated using Eq. (|6]l or O, 
depending of the direction of the forcing EMF. Its value is re- 
ported in column 9. 

Finally, we foUow iGuan & Gammid ( |2009|) and associate a 
turbulent "viscosity" v, with the flow, where 



(10) 



which we use to define a turbulent magnetic Prandtl number 
Pnit - Vt/rj, given in column 10. 

4.1. Radial diffusion of a vertical magnetic field 
4.1.1. £()v =4x 10"'° 

In this section, we describe in detail the results we obtained for 
models Ey-4E-10-R, Ey-4E-10-R-LONG and Ey-4E-10-A. All 
were obtained using a forcing EMF directed along the azimuthal 
direction, and varying only in the radial direction: Eq^ = and 
£■(),, - 4x10"'". We used Ae,^ - H. The first two were computed 
using Ramses and differ only in the duration of the averaging 
procedure we applied. The third was computed using Athena 
and serves both to validate our method and to quantify the un- 
certainty in our estimate of the turbulent resistivity and magnetic 
Prandtl number. 

The results we obtained for this set of parameters are illus- 
trated on figure|5]for model Ey-4E-10-R-LONG. There are four 
panels. The solid line in the first three (upper left, upper right and 
lower left) plots the radial profile of B ^, By and B. respectively. 



averaged in time between f = 30 and f = 150 using 120 dumps, 
and averaged in space over y and z. All plots use the same verti- 
cal scale. It is apparent from the lower left panel that the toroidal 
forcing EMF created a sinusoidally varying vertical field, as ex- 
pected from the discussion in section |2l The dashed line shows 
the sinusoidal curve obtained by a least square fit to the data, 
and gives an amplitude Bq- ~ 6 X 10"^ (this translates into a 
value j6,„av = 5.5 x lO'* as shown in Table |2l). It is immediately 
obvious than this rather low value is mostly due to the effect of 
turbulence. Indeed, using Eq. (|5]l, the effect of the microscopic 
resistivity alone would give an amplitude B^^ ~ 8 x 10""*, larger 
by about two orders of magnitude than the measured value of the 
vertical field. Using Eq. (|6]l, we converted Bo,- into a turbulent re- 
sistivity Tjiurb.x = 1.06 X 10"^. Combined with the time averaged 
value of ff = 2.0 x lO""*, this gives a turbulent magnetic Prandtl 
number of fm, = 1.27. 

Given the rather low value of the forcing EMF, it is reason- 
able to question the accuracy of these measurements. Indeed, the 
second and third panels in figure |5] indie ate that the fluctuations 
in By (upper right panel) are of comparable amplitude to Bofl 
The fourth (lower right) panel, which shows the time history of 
the rms fluctuations in the running time and spatial averages of 
B, (solid line) and B,, (dashed line), indicates a similar trend: al- 
though the rms fluctuations of B, become larger than those of B,. 
after about 20 orbits (thus indicating than the standard averaging 
duration of 50 orbits we usually take in the following is enough), 
the ratio between the two is only a factor of about two after 120 
orbits. In other words, the vertical field generated by the forc- 
ing EMF is not much larger than the turbulent fluctuations in the 
field. As such, our measurements are susceptible to a rather large 
uncertaintyQ In order to quantify this uncertainty, we consider 
model Ey-4E-10-R, which is the same as model Ey-4E-10-R- 
LONG but averaged over only 50 orbits, and model Ey-4E-10-A, 
identical to model Ey-4E-10-R but performed using Athena. As 



' The absence of fluctuations in the y and z averaged radial profile of 
B , comes from the solenoidal nature of the magnetic field 

- The reason to study such small fields is so that they do not signifi- 
cantly affect the property of the turbulence itself (the wavelength of the 
most unstable wavelength of the MRI associated with the field gener- 
ated by the forcing EMF is largely underesolved in our simulations). 
Thus, they merely act as a passive probe of the turbulence, which is not 
the case for the larger amplitude of E we consider in section |4. 1.21 
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Fig. 5. Time averaged radial profile of the vertically and azimuthally averaged radial (upper left panel), azimuthal {upper right 
panel) and vertical {lower left panel) magnetic field for model Ey-4E-10-R-LONG. On the third panel, the solid line shows the curve 
extracted from the numerical simulations while the dashed line is a least-squared fit to the data assuming a sinusoidal profile. The 
lower right panel shows the time variation of the standard deviation of the three magnetic field components, plotted using solid, 
dotted and dashed lines respectively. These plots illustrate that the forcing EMF creates, as expected, a radially varying vertical field 
with an amplitude larger than the fluctuations of the other magnetic field components 



shown on Table |2l we found extremely close values for the tur- 
bulent resistivities for all three models. The measured turbulent 
magnetic Prandtl number varies by about 50%, from Pm, ~ 1 to 
Pnit ~ 1.5. This diff'erence can immediately be attributed to vari- 
ations in the measured a, ranging from 1.7 x 10"^ to 2.3 x 10"^. 
Such a variation is compatible with the standard deviation of a 
reported in Table [1] and we conclude therefore that a safe esti- 
mate of the turbulent Prandtl number in the radial direction in 
this case lies in the range [1., 1.5]. 

This esti mate is in rough a g reeme nt wi th the recent results 
report ed by iGuan & Gammid (l2009l) and iLesur & Longarettil 
(|2009|) in simulations performed in the presence of a net mag- 
netic field in the azimuthal or vertical direction. It is also re- 
assuring that our measured value for the turbulent resistivity is 
in rough agreement with mean field theo ries like the Second- 
Order Correlation Approximation (SOCA. Ikadler & Rheinhard^ 
|2007|) also known as the First Order Smooth ing Approximation 
(FOSA, iBrandenb urg & Subramanianll2005h . In this approach, 
the turbulent resistivity is expected to take the form 

'7^^'''^=5v2w, (11) 

where Tcorr is the correlation time of the turbulent velocity fluc- 
tuations. Estimates for Tcorr are not straightforward to obtain. In 
a set of local and global str atified simulations of MR I-induced 
turbulence having a ~ 0.0 1, lFromang & Papaloizoul ( |2006|) and 
iFromang & NelsorJ (l2009t) found Tcorr ~ 0.157^0^;,. Adopting this 



estimate here, and using dv^ ~ O.llco given in table |2] for the 
three models we consider in this section, Eq. ( fTTT ) then gives 

T/^f'^'* = 1.14X 10"^ (12) 

a value that is, given the approximations involved in the SOCA, 
in very good agreement with the results reported in Table |2] 

4.1.2. Varying the forcing ElVIFs 

In the previous section, we considered a forcing whose ampli- 
tude only produces a very weak effect on the underlying turbu- 
lence. In this section, we relax this hypothesis by considering 
larger forcing EMFs. Such forcing may have an effect on the 
saturated properties of the turbulence itself, which in turn may 
change the values of the turbulent resistivity and Pm,. 

In the run we performed, listed in Table |2l the amplitude of 
the EMF was varied between l^l = 4 x 10"'" (described in the 
previous section) and l^l = 8 x 10"^, keeping Ae,.x = H in all 
cases. We tried larger values as well but found that, regardless 
of the code we used, the flow developed regions of low den- 
sity and large magnetic field. In these regions the large Alfven 
speed resulted in an extremely small timestep, which prevented 
the simulations from being continued. 

Figure |6] shows the radial profile (time averaged for 50 or- 
bits and spatially averaged over the azimuthal and vertical di- 
rections) of the vertical magnetic field we obtained in models 
Ey-lE-9-R {left panel), Ey-4E-9-R {middle panel) and Ey-8E- 
9-R {right panel). As was the case for the lower left panel of 
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Fig. 6. Vertically, azimuthally and time averaged radial profile of the vertical magnetic field for models Ey-lE-9-R {left panel), 
Ey-4E-9-R {middle panel) and Ey-8E-9-R {right panel). On each panel, the dashed line is a sinusoidal fit to the simulation data used 
to estimate the turbulent resistivities reported in Table |2] 



5? 2x10 




Athena (model Ey-4E-9-A in Table|2]i. The results are seen to be 
consistent with the results obtained with Ramses. 

Since Eq. ( fTTT i predicts the scaling of 77, j with the velocity 
fluctuations, in figure [T] we show our results in the {Sv^, 7], ^) 
plane. The data obtained with Ramses are plotted using dia- 
monds while the points computed using Athena appear as a star 
symbol. The dotted line is a naive estimate that uses the pre- 
diction of the SOCA given by Eq. ( fTTT i. assuming that the cor- 
relation timescale is independent of the velocity fluctuations: 
Tcorr - Q-^iTgrh- Clearly, the dotted line does not fit the data. 
This is most likely because Tcorr varies as the strength of the tur- 
bulence changes. On dimensional grounds, it can be estimated 
that 



^x,corrl ■> 



(13) 



Fig. 7. Turbulent resistivity in the radial direction as a function 
of the turbulent velocity fluctuations for various models com- 
puted using different amplitude of the forcing EMFs. Diamonds 
correspond to run performed using Ramses while stars use the 
results of models computed with Athena. The dotted line shows 
the prediction of SOCA assuming a constant correlation time, 
independent of the velocity fluctuations, while the dashed line 
shows the same prediction but computed assuming a correlation 
lengthscale for the magnetic field independent of the amplitude 
of the velocity fluctuations {see text for details). The latter is seen 
to be in good agreement with the data. 



figure |5] the solid line plots the results of the simulation and the 
dashed line shows the sinusoidal fit used to derive the turbulent 
resistivity and magnetic Prandtl number. As expected, figure |6] 
shows that the amplitude of the magnetic field increases with 
the forcing amplitude. It also illustrates the quality of the sinu- 
soidal fit for these models. The turbulent resistivity and magnetic 
Prandtl numbers we calculate using the amplitude of the steady- 
state field are shown in Table |2] We also report (columns 4 to 
7) the time averaged a values and time averaged velocity fluctu- 
ations measured in these runs. As forcing is increased, we find 
that Tjt x increases. Increasing the forcing by a a factor of 20 in- 
creases rjix by about 2.5. At the same time, a roughly triples. 
This leads to a mild but systematic increase in Pm, with the am- 
plitude of the forced vertical magnetic field. To further test our 
results, we also computed a model with l^l = 4 x 10"^ with 



where lx,corr is a correlation length associated with the verti- 
cal magnetic field fluctuations. Although it is rather ill-defined, 
an order of magnitude es timate for lx,corr can be obtained as 
(iLesur & Longarettil2007h 



j B,{x,y,z)B^{x' ,y,z)dxdx' 
J Bl{x,y,z)dx 



(14) 



where angled brackets denote a spatial average over the y and z 
directions followed by a time average. Using Eq. (fT4l i. we found 
that lx,corr is Only weakly varying with the strength of the turbu- 
lence. Indeed, we obtained 0.06// < h.corr < 0.07// for forcing 
amplitudes ranging from l^l = (i.e. no forcing) to l^l = 8x10"' 
(i.e. the largest forcing amplitude we considered). It would be 
dangerous at this stage to take the estimate given by Eq. ( fT4l ) as 
a good way of calculating Tcorr, as Eq. ( fT3] ) is only valid on di- 
mensional grounds. Instead, we used the independence of lx,corr 
to write the correlation timescale as 



5v 



To- 



6vx 



(15) 



where tq - Q.\5To,-i, and Svxq denotes the velocity fluctuations 
in the absence of forcing (given in Tabled}. This gives 



SOCA c c 



(16) 



which we used to compute the dashed line shown in figure|2] The 
agreement with the data is much better and captures the scaling 
of the turbulent resistivity as the amplitude of the forcing EMFs 
is varied by more than one order of magnitude. 
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4.2. Radial diffusion of a vertical field in extended boxes 

Next, we consider the turbulent resistivity in boxes extended in 
the radial direction. These models are based on the results of 
model A-LB, which we restarted at time t - 30 with forcing of 
different amplitudes. We only consider forcing in the azimuthal 
direction in this case. Again, we report the properties of the tur- 
bulence and the values of 77, j and Pm, we obtained in Table |2l 

To check the results obtained in smaller boxes, we first con- 
sider model Ey-4E-9-A-L4-lxl . In this case, we used Ae j^ - H 
and EQy = 4 x 10 Thus, it is the big box equivalent of mod- 
els Ey-4E-9-R and Ey-4E-9-A. We therefore expect it to yield 
the same results. Table|2]shows that this is indeed the case, with 
the values of a, rj,^ ^ and Pm, being nearly the same in all three 
models. 

Next, we turn to models Ey-4E-10-A-L4-lx4 and Ey-4E-9-A- 
L4-lx4 for which Ae,x - 4//. The first result to note is that the 
amplitude of the vertical field we obtained is larger in these two 
cases than in the small box models with equal forcing EMFs. 
This was somewhat expected. Since Ae,x is now larger, the lo- 
cal gradient of the vertical field are smaller in the big boxes and 
the diffusion of magnetic field is less efficient. Nevertheless, the 
measured turbulent resistivity are significantly larger than their 
small boxes counterpart. At the same time, the velocity fluctua- 
tions in the radial direction are comparable to the small boxes 
models having the same forced vertical field amplitude. It is 
likely that all other statistical properties of the turbulence share 
the same insensitivity to box size, which means that a straight- 
forward application of the SOCA breaks down in this case. The 
increased resistivity is therefore due to the increased wavelength 
^E,x of the forcing EMF used in that case. One possible explana- 
tion, to be taken with care, for such large 77, can then be described 
as follows: the resulting forced magnetic field, having a larger 
scale, will most hkely be diffused by larger eddies. Because the 
kinetic energy power spectrum is decreasing with wavenumber, 
such eddies have a larger kinetic energy that makes them more 
efficient at diffusing the magnetic field. 

The consequence of the larger turbulent resistivities obtained 
in our large box simulations is lower values of Pm,. In the case of 
model Ey-4E-10-A-L4-lx4, Pm, ~ 0.5, i.e. a factor of about three 
times smaller compared to the small boxes. We recover Pm, of 
order unity for model Ey-4E-9-A-L4-lx4 owing to a larger value 
of a. This is presumably due to the fact that the flow starts to be- 
have as if threaded by a net vertical flux for this long wavelength, 
large amplitude forcing. 



4.3. Vertical diffusion of an azimuthal magnetic field 

In addition to studying radial diffusion of vertical magnetic fields 
as described in the previous section, we also performed simula- 
tions in which we study the vertical diffusion of an azimuthal 
magnetic field. In this case, we restricted our analysis to small 
boxes, and we chose Eoy - and therefore considered forcing 
EMFs purely in the radial direction. We allowed for variations 
of ffie amplitude of the forcing in the range 4 X 10"' < Eq^ < 
8 X lO"**. We also tried Eq^ - 4 x 10"'" as in the case of radial 
diffusion. However, because the azimuthal magnetic field fluc- 
tuations are of larger amplitudes than the vertical field fluctua- 
tions, it has proven difficult to extract a meaningful mean field 
from that model. We therefore decided to consider only large 
amplitudes for the forcing EMFs. 

The results we obtained are summarized once more in 
Table|2] We obtain turbulent resisitivities that are generally lower 
than for the case of radial field diffusion. This is not inconsistent 



with the prediction of SOCA, as the velocity fluctuations in the 
vertical direction are of lower amplitudes than the velocity fluc- 
tuations in the radial direction. Using 



SOCA 
lt,z 



(17) 



we obtained tjI'^^'^ = 6.8 x lO"*" for model Ex-4E-9-R, while 
Table |2] gives 7;,^ - 1.05 x 10"^ (as compared to ri, x = 1.81 x 
10"^ for model Ey-4E-9-R). Although it is less accurate than the 
estimates we obtained in the case of radial field diffusion, the 
prediction of the SOCA still gives a reasonable value for the 
turbulent resistivity. 

Finally, because of the lower values of the turbulent resis- 
tivity, the final column of Table |2] reports systematically larger 
values of the turbulent magnetic Prandtl number, ranging from 
1.57 for model Ex-4E-9-R to 3.00 for model Ex-4E-9-R. As in 
the case of a vertical field diffusing radially, Pm, is a slowly in- 
creasing function of the strength of the turbulence. 



5. Discussion and conclusion 

In this paper, we have presented a set of local numerical simu- 
lations aimed at measuring the properties of turbulent diffusion 
of magnetic field, quantified using an anomalous turbulent resis- 
tivity 77,, resulting from MHD turbulence induced by the MRI. 
We considered only the case in which turbulence develops in 
the absence of a mean magnetic field. We used two different 
codes, Athena and Ramses, and vary both the box size and the 
magnetic field geometry. The main result that emerges from our 
simulations is that 77, tends to be slightly smaller, but similar 
in magnitude, to the anomalous viscosity v, that can be asso- 
ciated with the outward transport of angular momentum due to 
the turbulence. In other words, their ratio, the turbulent magnetic 
number Pm,, is of order unity (actually in most of our runs it 
was slightly larger than one). We also found that our results are 
roughly consistent with the predictions of the mean field theory 
SOCA. Perhaps more relevant is the fact that our results are con- 
sistent with those recen tly published by Guan & Gammie (200^ 
and lLesur & Longarettil (120091) . who used different field topolo- 
gies (toroidal and vertical mean field), different numerical meth- 
ods (similar to those in the ZEUS code, and an incompressible 
pseudo-spectral code respectively) and, as described in the in- 
troduction, different approaches to measure the turbulent resis- 
tivity. This broad agreement demonstrates that all these results, 
obtained and published independently, are largely insensitive to 
the numerical method and the topo goly of t he magnetic field. As 
suggested already in the literature (Parker 1971), a safe first or- 
der guess of the magnitude of the turbulent resistivity is therefore 

77, ~ V,. 

A number of other results emerge from the comparison 
of the suite of simulati ons we present here wi th the work of 
iGuan & Gammiel(l2009h and lLesur & L ongaretti (2009). First, it 
is wor th noting that our results and those of Les ur & Longarettil 
(l2009l) were obtained inclu ding explicit dissipation, while those 
of iGuan & Gammi3 (1200 9^ relied only on numerical dissipation. 
Since all the results are roughly consistent with one another, this 
is apparen tly not an issue as far as dete rmining Pm, is concerned. 
Moreover. iLesur & Longarettil (l2009h consider the case Pm - 1 
while we have used Pm - 4. The broad agreement between both 
studies therefore suggest that Pm, does not strongly depend on 
microscopic dissipation coefficients. However, it would be pre- 
mature to draw definite statements here since other differences 
between both works might mask a possible Pm dependence in 
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the results. Although very computationally expensive, a system- 
atic study of the sensitivity of our results to the value of Pm 
is needed in the future. Second, Pm, appears to be smaller in 
large boxes. Indeed, we found Pm, = 0.57 for model Ey-4E- 
10-A-L4-lx4, while Pm, is found to be larger than one in small 
boxes having identical parameters. This rather low value, ob- 
tained in a large box and for a small forcing EMF, is very close 
to the values quoted by Guan&i Gamrnie (2009) for models hav- 
ing t he same radial box size a nd a weak imposed field. In addi- 
tion, iGuan & Gammid (12009 *) considered boxes larger than 4H 
in the radial direction and found this value of Pm, to be inde- 
pendent of box size . By c ontrast, both our results and those of 
iLesur & Lonparettil (l2009l) suggest a value larger than one (or 
equivalently fairly lo w values for the t u rbulen t resisitivity) in 
small boxes, a regime iGuan & Gammid (l2009l) did not investi- 
gate. Taken together, all these results point toward a decrease of 
Pm, when box size is increased from values of a few to values of 
roughly one half, with a plateau obtained for boxes larger than 
4H. 

Despite their broad agreement, it is instructive to consider 
the di fferences between our results and the aformentionn ed pa- 
pers dGuan & Gammid 120091: iLesur & Longarettil l2009t) . One 
such difference is that, regardless of the model or the box size, 
we alway s find that Pm, increases wi th the amplitude of the forc- 
ing EMF. iLesur & Longaretti (2009) also report such a trend for 
the models they label XXZn, in which they consider the vertical 
diffusion of a toroidal field, as we do in section l431 These trend 
is not so clear for their other cases. However, most of their runs 
were obtained in the presence of a net vertical flux. The nature of 
the turbulence is strongly modified in that case and this is prob- 
ably not appropriate to d raw the comparison with those models 
any further. By contrast, iGuan & Gammid (l2009h report no de- 
pendence of Pm, on the amplitude of the magnetic field they im- 
pose for models that are similar to ours (i.e. the diffusion of a ra- 
dially varying vertical field). Even if the different field topology 
they consider (namely a net toroidal flux) plays a role, the differ- 
ence with our results can most likely be attributed to the differ- 
ent methods that are being used. While the present paper aims to 
measure a steady state response of the flow to a_perturbation that 
is imposed at all times, iGuan & Gammid (|2009j) Superpose an 
additional field to an already turbulent flow at f = and let it de- 
cay. Thus it is possible that transients associated with this addi- 
tional field may complicate the estimate of time averaged values 
for the transport coefficient^ In particular, for large strength of 
the additional field, the flow may not have enough time to reach 
fully saturated values of the turbulence before the strength of the 
imposed field decays significantly. This might lead to an under- 
estimate of the value of a and consequently of the value of Pm,, 
masking the increase of Pm, with forcing that we found. 

It is also important to stress that all three approaches share 
common limitations to their analyses. First and foremost, the 
analysis is by definition local, whereas the diffusion of mag- 
netic field accross the disk is an intrisically global problem 
that depends on the large scale properties of the disk such 
as the radial profiles of th e surface density and magnetic flux 
dSpruit & Uzdenskvll2005h . Another related limitation is the ne- 
glect of vertical density stratification in the disk. Indeed, all three 
studies presents numerical simulations performed in "unstrati- 
fied" shearing boxes, neglecting the vertical component of grav- 



^ Of course, such transients are most like ly of physical origin, in the 
sense that the decay timescale measured bv lGuan & Gamniig l l2009h is 
related to the relaxation timescale of the turbulence. Nevertheless, their 
presence complicate the estimate of Pm, 



ity. In simulations including gas density stratification, it was 
found that MHD turbulence is suppressed in the lo w density 
corona above and below the disk midplane dStone e t al. 199^. 
It is likely that turbulent diffusion will be greatly reduced at 
those locations, possibly causing the magnetic flux to remain 
anchored to the disk corona. This could prevent diffusion that 
otherwise would have been driven by turbulence in the mid- 
plane, or by channel modes that couple the upper layers of 
the disk with the midpl ane (Bisnovatvi-Kogan & Lovelacel2007l: 
iRothstein & Lovelacdl2008; .Lovelace et al...2009l) . A numerical 
test of this scenario is beyong the scope of the present paper but 
should be consider ed in the future. Finally, let us mention the 
very recent work of iB eckwith et al.l (12009 ). Using global simu- 
lations that combine the large scale nature of magnetic field dif- 
fusion and the vertical density stratication (but, of course, with a 
dramatic reduction in the spatial resolution), these authors chal- 
lenge the concept of turbulent resistivity itself and advocate a 
completely different scenario to describe magnetic flux evolu- 
tion in accretion disks. Clearly, despite the agreement with al- 
ready published calculations, the results presented here should 
be taken with care when used in the context of astrophysical ap- 
plications, such as jet launching or magnetic flux distribution in 
accretion disks. 
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